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Abstract 



Focussing on preplanetary grains growth, we discuss the properties 
of dust aggregation driven by magnetic dipole forces. While there is no 
direct evidence for the existence of magnetic grains present in the solar 
nebula, there are reasons to assume they may have been present. We 
derive analytical expressions for the cross-section of two interacting 
dipoles. The effective cross section depends upon the strength of the 
magnetic dipoles and the initial velocities. For typical conditions the 
magnetic cross section is between 2 and 3 orders of magnitude larger 
than the geometric cross section. We study the growth dynamics 
of magnetic grains and find that the mass of the aggregates should 
increase with time like t^'^ whereas Brownian motion growth behaves 
like t^. A numerical tool is introduced which can be used to model 
dust aggregation in great detail, including the treatment of contact 
forces, aggregate restructuring processes and long-range forces. This 
tool is used to simulate collisions between magnetic grains or clusters 
and to validate the analytical cross-sections. The numerically derived 
cross section is in excellent agreement with the analytical expression. 
The numerical tool is also used to demonstrate that structural changes 
in the aggregates during collisions can be significant. 



1 Introduction 

Coagulation of dust grains is generally believed to be the mechanism by 
which growth of particles in the early solar nebula proceeds from sub-micron 
grains all the way up to planetesimals (|Weidenschilling and Cuzzi, 1993| ; 



Beckwith et al., 20001) . However, many questions remain about the detailed 



way in which coagulation proceeds, in particular during the initial phase when 
the particles are still very small and dynamically well coupled to the gas. If 
the particles are very strongly coupled to the gas, collisions have to rely on the 
relative velocities induced by Brownian motion ( [Kempf et ai, 1999| ). While 
these slow collisions are certain to result in almost perfect sticking between 
the grains QChokshi et ai, 1993| ; |Poppe et ai, 2000| ), it turns out that they 



may be too slow to produce the necessary dust growth in the limited time 
frame available ( |Kempf et ai, 1999|) . First of all, the low relative velocities 



severely limit the number of collisions taking place. As a second complication, 
the energies involved in the collisions are insufficient to change the structure 
of the aggregates (|Dqminik and Tielens, 199"^ ; [Blum and Wurm, 2000|) . 
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Simple hit-and-stick growth in a cluster-cluster aggregation process leads 
to extremely fluffy aggregates with low fractal dimensions ( [Meakin, 1991 



Kempf et al., 1999| ). |Kempf et al.\ showed that the Brownian motion 



leads to fractal dimensions D < 2 for which significant changes in the cross- 
section-to-mass ratio can only be expected for very large particles. The 
friction time, i.e. the time in which a particle adapts its velocity to the 
surrounding gas, changes in a Brownian growth process only as r^/^ with 



the particle radius. Furthermore, [Kempf et al.\ neglected the rotation of 



the aggregating particles which will lead to even lower fractal dimensions of 
the collision product ( [Blum et al., 20"00| ). Since the aerodynamic behavior 
of particles is governed by the friction time ( Weidenschilling, 1977 ), the 
aerodynamic properties of large fluffy aggregates and of small grains are 
almost identical. This poses a serious problem for the coagulation process 
which must accelerate eventually, in order to fulfill the time constraints given. 

Possible ways to accelerate growth have been discussed in the literature, 
in particular grain settling and radial drift (e.g. Weidenschilling, 1980 ), or 
differential coupling to turbulent eddys of different sizes ( [Weidenschilling j 
1984[ ; [MizuNO et al, 1988] [Markiewicz et al, 1991[) . Furthermore, aerody- 



namic concentration is considered as a mechanism to increase densities locally 
for enhanced coagulation ([Cuzzi et al, 2001[ ; [Klahr and Henning, 1997[) . 
However, most of these processes rely on a diversity in the aerodynamic 
properties of particles: Settling or radially drifting grains will only induce 
enhanced collision rates if not all grains settle or drift at the same rates. 
Turbulent eddies will only enhance collision rates if different grains couple to 
different eddie scales. [Weidenschilling (1997| ) reports that starting from 
10~^cm grains, a hight central density of settled dust can still be achieved 
even if the fractal structure of the grains is taken into account. 

Coagulation may proceed on quite different paths if the basic mechanism 
of coagulation is different from what is normally assumed. An interesting idea 
was put forward by [Nuth et al. (1994[ ) and [Nuth and Wilkinson (19951 ) 
who studied aggregation of magnetic grains in laboratory experiments. They 
found that this accretion process proceeds rapidly forming large networks 
of linear aggregates. Nuth et al.\ suggested that such networks may be 
important for grain aggregation in the solar nebula. One way in which this 
mechanism could change the course of aggregation is to turn a cluster-cluster 
aggregation into a particle-cluster aggregation process. If the networks of 
magnetic grains can form quickly enough, they may eventually provide a 
large fraction of the available coUisional cross-section in the nebula. The main 
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aggregation process for non-magnetic grains would then be the colhsion of a 
grain with such a network. This mechanism would be physically similar to a 
cluster-particle aggregation process which produces aggregates with fractal 
dimension Dj = 3 ( [Ball and Witten, 198^ ), leading to more compact 
structures quickly. 

In this article, we will present the foundations for a comprehensive investi- 
gation of magnetic dust aggregation from a theoretical and an experimental 
point of view. We start by reviewing the evidence for magnetic grains in 
the solar nebula. We will then focus on the physical implications of dipo- 
lar interaction between individual dust grains, i.e. enhanced collisional cross 
sections and accelerated aggregation dynamics. We also present a new nu- 
merical tool we devised for the investigation of dust aggregation including 
a detailed treatment of grain-grain interaction through mechanical contacts, 
and dipolar (magnetic) forces. We carried out numerical simulations of 2- 
body collisions with magnetic grains in order to test the code and compare 
the results with the analytically derived cross sections. Two subsequent pa- 
pers will communicate our experimental and numerical work of magnetic 
grain growth on the basis of the material presented here. 



2 Magnetic dust in the solar nebula 

If magnetic interaction is to change the dust aggregation scenario, we have 
to assume the existence of a nebular dust component carrying remanent 
magnetization (|Nuth et al, 1994| ; [Nubqld and Glassmeier, 1999| , |2000|) . 



The following paragraphs present evidence as to the extent this assumption 
seems to be justified. 

Formation of chondrules Many of the small number of pieces that have 
been found to complete the puzzle of solar nebular physics come from the 
study of meteorites. [Kerridge and Matthews (1988| ) and [Kerridge 
|(1993|) provide a comprehensive overview of this subject. 

An important result from the study of meteoritic material is the occur- 
rence of chondrules, i.e. small spherical inclusions which are embedded in 
a presumably more primitive matrix. Despite a large number of different 
hypotheses ( [Boss, 1996| ) as to the origin of chondrules, the consequences 
of the necessary melting process have been studied intensively. Chondrules 
are formed when a physical process in the solar nebula melts silicate dust 
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aggregates followed by rapid cooling of the molten material which congeals 
as spherical droplets ( |Gibbard et al, 1997|) . Nebula lightning, shock fronts, 
and friction have been discussed as possible energy sources. The heating 
process itself has been repeatedly investigated both theoretically and experi- 
mentally (ILevy, 19881 ; [Grossman et al, 198^ ). The heated dust aggregates 



of possibly interstellar composition undergo internal reduction-oxidation re- 
actions, which result in an evaporative loss of metallic iron from the melt 
(|NUTH, 1989|) . The iron vapor then condenses into sub-micron- sized dust 



grains. Depending on the grain size these grains can have different magnetic 
properties ( |J\uth et al., 1994| ; [Withey AND JNUTH, 1999| ) . Very small grains 



below 20 nm are superparamagnetic, large grains > 80 nm show multidomain 
behavior. Intermediate sizes, in particular grains between 40 and 60 nm are 
single domain grains which spontaneously become magnetic dipoles. In a 
more general way, melting events in the solar nebula will always produce 
magnetic minerals such as magnetite ( [Rietmeijer et al, 1999|) . 

Chondrules themselves originate from already existing aggregates, so the 
iron grains produced by this process will not have been present very early on 
in order to start the aggregation process. However, similar melting events 
must have been a frequent phenomenon in the early solar system. The as- 
sumption of an important population of magnetic iron dust appears plausible 
on these grounds. 



The paleomagnetic record in meteorites As summarized by [SUGIURA 



AND Strangway (1988D , many meteorites exhibit natural remanent magne- 
tization (NRM). The carriers of the NRM are magnetic minerals and metals 
such as iron-nickel (Fe-Ni: kamacite, taenite), phyrrotite (Fei_xS), and mag- 
netite (Fe304; Jackel, personal communication). A possible source for the 
NRM in meteorites are strong magnetic fields that are assumed to have ex- 
isted in the solar nebula, magnetizing the dust as it cooled through the Curie 
point. 

It is usually assumed that the NMR in meteoritic material was accuired 
after the formation of the meteorite parent body. However, some meteorites 
show evidence of heterogeneous magnetization, which could be a hint of mag- 
netization that existed before the dust was incorporated into the meteorite 
parent body. The magnetization of isolated grains is certainly more easily 
explained than the magnetization of macroscopic objects, for the latter would 
require very stable geometric conditions, i.e., a fixed orientation of the ob- 
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ject that is to be magnetized relative to the magnetizing field. In any case, 
the paleomagnetic record in meteorites also points toward the existence of a 
magnetized dust population in the solar nebula. 

Cometary dust Cometary nuclei are widely supposed to contain the orig- 
inal building material of our solar system in almost pristine form. The 
great effort which is put into prestigious and adventurous space missions 
like Rosetta, Stardust, and Deep Space 4/Champollion is a measure of the 
importance attributed to the analysis of cometary matter. 

During the Giotto and Vega missions to comet p/Halley in 1986, mass 
spectrography of cometary dust grains impacting the spacecraft gave some 
insight as to the chemical composition of cometary dust. [SCHULZE et al\ 
(1997|) summarized the results of these measurements. Accordingly, Halley's 
dust is composed mainly of magnesium-rich silicates, but there is an impor- 
tant fraction (> 10%) of potentially magnetic minerals and metals similar 
to those which have been identified in meteorites, i.e. kamacite (nickel-iron), 
phyrrotite (Fei_xS), and magnetite (Fe304). 

Some interplanetary dust particles (IDPs) appear to originate from comets 
( |Brownlee, 1985]) . Up to now, no measurements of the magnetic properties 
of IDPs have been obtained (Bradley, priv. comm.). The presence or absence 
of aggregates of magnetic grains in cometary IDF would be a direct way to 
determine if a large population of magnetic grains in the outer solar nebula. 

Asteroids On its way to the Jovian system, the Galileo spacecraft inves- 
tigated the asteroids Ida and Gaspra. The magnetic field measurements in 
the vicinity of the two bodies show significant changes in the solar wind 
magnetic field as it encounters the asteroids. The magnetic signatures have 
been interpreted in terms of a remanent magnetization of the two asteroids 
by [KlVELSQN et al. (1993| ). The same holds for the encounter of the Deep 
Space 1 (DSl) spacecraft with the asteroid Braille in 1999 ( |Richter, 2001|) . 

Since asteroids - which should have sizes similar to planetesimals during 
solar system formation - are much too small to sustain a dynamo, their 
magnetization could be the trace of primordial magnetic fields or the result 
of the accretion of previously magnetized material. 

Dust around young stars ISO observations of young nearby stars have 
been used to derive the dust properties in the disk around stars close to 
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the main sequence. Unfortunately, most T Tauri stars are too faint for this, 
but their higher-mass counterparts, the Herbig AeBe stars, have proved to 
be excellent targets for such a study [Waters and Waelkens (1998| ). In 
several of these stars, the signatures of crystalline silicates have been found 
(e.g. IMalfait et a/., 19"98| ; |Bouwman et a/., 2000|) which are very likely 



formed by thermal annealing ( [Bouwman et al., 2001| ). These silicates are 
magnesium-rich, which may leave the iron free to form metallic iron and 
other potentially magnetic grains. 

3 The influence of magnetic interaction on 
the dust aggregation process 

In this article we focus on magnetic forces as the only long-range force. Elec- 
trostatic interactions between charged grains could be another source of such 
long-range interactions between dust particles. Their effect on dust coagula- 
tion in a nebular context has been investigated both theoretically ([Horany] 



AND GOERTZ, 1990| ; pssENKOPF, 1993| ) and experimentally ( [Marshall 



[and Cuzzi, 2001| ). Since the solar nebula is generally considered as opti- 
cally thick at least in its central regions, photoelectric charging can be ruled 
out when discussing possible charging mechanisms for dust grains. Triboelec- 
tric charging would be due to collisions which usually lead to direct sticking 
within the range of relative velocities considered here, even if van der Waals 
interactions are the only attracting force. Therefore, electrostatic interac- 
tions are probably of importance only for larger grains ( Marshall and 
Cuzzi, 20011) . 



In the absence of other long-range forces, random thermal agitation is the 
sole competitor for magnetic interaction. The relative importance of these 
two effects is generally given by a dimensionless parameter A defined by: 

AMbT' ^ ' 

where /i is the magnetic moment of the particles, /iq is the magnetic field 
constant and fc^ is the Boltzman constant. The parameter d in Eq. ^ defines 
a certain characteristic interaction distance which is usually taken to be equal 
to the particle diameter, d = 2R, or to the sum of the particle radii, d = 
Ri + Rj, if we are dealing with multi-disperse grains. For single-domain 
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particles, which are typically micron-sized or smaller, A can take on values 
of the order of 10^ at T = 100 K. In these cases, the aggregation process will 
not be diffusive but is dominated by magnetic forces. Furthermore, magnetic 
aggregates should be stable against destructive thermal influences. 

According to statistical mechanics, dust particles should acquire a ther- 
mal rotational energy f/jot = which amounts to very high angular 
velocities for micron-sized particles. In order to have magnetic interaction 
overcome this thermal spin-up, the particle-particle interaction energy must 
be at least of the same order of magnitude. [Nubold and Glassmeier] 
(2000| ) have shown that this happens inside an interaction radius given by: 

= 2d\"\ 

which is about 10 to 100 times larger than the actual particle radius for 
strongly magnetized grains. It is on these length scales that magnetized dust 
particles "see" each other and interact. In the following sections we will 
examine some of the consequences in more detail. 

3.1 Collisional cross section 

Grain-grain collisions and sticking in the dilute nebular gas result from rel- 
ative motion, e.g. due to a thermal distribution of grain velocities. Owing 
to the minute geometric dimensions of the dust particles, the collision prob- 
ability is usually very low. Long-ranged forces, e.g. magnetic forces between 
magnetized grains, lead to an enhancement of collisional cross sections. 

Usually, the collisional cross section of a (spherical) particle is simply 
given by its geometric projection area, i.e. 

a = nR^, (3) 

where R denotes the particle radius. For two colliding particles we have: 

a = 4 ■ 7i{R, + Rjf . (4) 

The collision probability in dust aggregation processes depends directly on 
the particle cross section: the bigger cr, the shorter the mean free time r ~ 
1/a between subsequent particle collisions. 
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NUBOLD AND Glassmeier (2000) investigated this effect for nanometer- 
sized iron grains numerically, based on the experimental results obtained by 
NUTH et at. (199^ . From scattering simulations with spherical iron particles 



with a diameter of 20 nanometers moving at thermal velocity, ^/k^TJm, it 
was shown that: 

— (^oo = t^th,)~ 900. (5) 
a 

Since the collision frequency depends on the cross section, an enhancement 
as given by Eq. |^ will cut aggregation time scales by almost three orders of 
magnitude. From Eq. ^ we see that the cross section may be of the same 
order as the magnetic interaction radius. 

In the calculations leading to Eq. |^ we implicitly assume that magnetized 
grains are always aligned with the ambient magnetic field, e.g. an external 
field or the field created by other magnetized dust particles. While this 
should be the case in terrestrial experiments due to high surrounding gas 
densities, it is not obvious that instantaneous relaxation of dipoles should 
occur in very dilute media such as the solar nebula. We will now relax this 
assumption. 

When a magnetic dipole ^ interacts with another magnetic dipole /x' over 
a distance r, according to [Jackson (1975| ) the interaction energy U can be 
calculated from 
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(6) 



The most important contributions to the interaction energy come from dis- 
tances r d and relative orientations 9ij 0, where Oij denotes the angle 
between two dipoles contained in a common plane. Thus, magnetic dipoles 
tend to align in a pole-to-pole configuration while minimizing the interaction 
energy: 

where we have introduced the potential parameter a. An analytical cal- 
culation - such as the Rutherford formula for electrostatic interaction - of 
the coUisional cross section between magnetized dipolar particles is not 
possible since we are not dealing with a central force problem where the 
interaction solely depends on the inter-particle distance r. In general, the 
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force between two magnetized dust grains will not be directed along the line 
connecting the particles centers. This is equivalent to saying that there is no 
conservation of angular momentum. However, Eq. |^ can be used to derive 
an approximate expression for a since it essentially describes a central po- 
tential, U = U{r), if we replace the variable d by r. Particle motion in this 
type of potential will be confined to a plane, and orbital angular momentum 
is conserved. The magnetic moments are assumed to be parallel to a line 
connecting the particle centers. 

Following the treatment of particle motion in a central potential (e.g. 
Landau and Lifschitz, 1984| ), scattering of two particles with mass m 



can be described in terms of an effective potential 

UMr) = U{r) + ^ , (8) 

in which we have introduced the reduced particle mass, m = m? /2m = m/2 
for identical grains. Equation ^ combines the interaction potential f/(r) given 
by Eq. |^ and the conservation of orbital angular momentum /. We are now 
dealing with an equivalent one-body problem. 

For scattering problems, I is usually given in terms of the impact param- 
eter h: 

I = m ■ ■ b, (9) 

where Voo denotes the relative velocity prior to interaction, e.g. thermal 
motion. Combining Eqs. |^ through |^ yields: 

(10) 

We refer to Eq. |^ for a definition of a. For potentials with r~^-dependence, 
the maximum impact parameter can then be derived by simply setting the 
initial energy 1/2 mf^ equal to the Ues as given by Eq. |r^. It should by 
noted, however, that the effective interaction potential for dipolar interaction 
(cf. Eq. |To|) does not display a local minimum as is the case for gravitational 
or electrostatic interaction, but a local maximum given by: 
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We therefore obtain femax by letting E = Uq, which imphes that an incoming 
particle first has to cross a potential wall of height Uq before it can fall into 
the center of force: 



2a2 V 3 J 2 
Therefore we find for the maximum impact parameter 



(12) 



3-(^)'''- (13) 



and the magnetic cross section is 



^m = 37r- ( ^] ' . (14) 



The central-potential approximation necessarily overestimates the magnetic 
cross section by a factor or order 2, because only maximum interaction is 
taken into account when deriving Eq. |1^. 

3.2 Accretional remanence 

A direct consequence of the increased particle cross section for magnetic 
grains is the selective coagulation of magnetized dust particles with their kin. 
Undisturbed magnetic aggregation should lead to the formation of chain- 
like configurations with very low fractal dimension, since magnetic grains 
tend to coagulate in a pole-to-pole pattern. Df = 1 characterizes a per- 
fect chain or rod. With non-magnetic dust grains, aggregation processes 
in three-dimensional Euclidian space usually result in structures with frac- 
tal dimensions around 2 ( [Kempf, 1997| ; |Wurm AND BLUM, 1998| ; |WURM 



19971 ; IWeidenschilling et ai, 1988|) . More recently however, [Blum et'o! 



(2000| ) and [BLUM AND WURM (20001 ) reported on the growth of aggregates of 



Si02-spheres with elongated structures and low fractal dimension (Df = 1.3 
(IBlum et ai, 20001) ) which can be attributed to rotational effects. 



In the case of magnetic grains, chain-like aggregation should coincide with 
an efficient conservation of the total magnetic moment. For perfectly aligned 
dipolar particles, the magnetic moment of a chain is given by the sum over 
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the individual magnetic moments /ij = /i: 

N 

/ichain = Ati = N ■ fi . (15) 

For sufficiently high values of the magnetic interaction parameter A (cf. 
Eq. |T|), i.e. A ^ 1, there will be no statistical mixing of individual mag- 
netic moments, but - leaving aside thermal agitation and locally disturbed 
agglomeration due to multi-particle interaction - an increase of the magnetic 
moment per cluster with the number of constituent particles. This increase 
Z^agg = f^agg{N) is called accretional remanence ( [Nubold and Glassmeier, 



1999|, pOOOp. Its functional dependence was determined from numerical sim- 



ulations ( [NuBQLD AND Glassmeier, 2000| ) to be 

fx,,,ifx,N)=fx-N'-'', (16) 

where 11 = Hi again denotes the magnetic (dipole) moment of a single grain. 
Large aggregates assembled from thousands of small magnetized grains will 
therefore still show a large fraction of the total magnetic moment of their 
constituents. 

[Figure 1 about here.] 
3.3 Aggregation Dynamics 

Due to its irreversible character, studies of the the kinetic properties of par- 
ticle aggregation focus on the dynamic properties of the distribution of clus- 
ters. Within this context, the features of aggregation are described in terms 



of the cluster-size distribution na{t) as a function of time (e.g. [Miguel and 



Pastor- Satorras, 1999|, and references therein). As is the case for all dis- 



tribution functions, Ug can be used to calculate higher moments, in particular 
the mean cluster size S{t) which is defined by: 

E^V(t) 

Sit) = — , (17) 

s 

the variable s being one particular cluster size from the distribution. The 
asymptotic behavior of the mean cluster size is a power law 

S{t) ~ e (18) 
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with the so-called dynamic exponent, z, which depends on the aggregation 
process. For ensembles composed of identical particles, the mean cluster size 
can be identified with the mean aggregate mass m{t). 

To our knowledge, the dynamical evolution of magnetic aggregation has 
never been analyzed before in an astrophysical context, i.e. involving micron- 
sized particles, low number densities, and a regime of free molecular flowQ. 
According to Plum et al. (20001) , we expect a behavior (cf. Eq. 



m{t) oc (19) 

for a quasi-monodisperse system. For an ensemble of aggregating particles, 
the dynamic exponent z depends on the evolution of the collisional cross 
section a oc and of the collision velocity Vc oc . The exponent z can 
be derived using the following equation ( |Blum et al., 200U| ; [Meakin, 1991] ): 



z — 1 , , 

v + jy = . (20) 

Just like m(t), Oc and Vc are assumed to obey power laws. Common dust 
aggregates have Df = 2, therefore dc oc m and rj = 1. Thermally induced 
relative velocities scale with m~^/^, i.e. v = —1/2. For such systems, we 
expect to find a dynamic exponent z = 2. 

In magnetically driven coagulation in free molecular fiow, particle motion 
is not diffusive but governed by strong inter-particle forces. For aligned grains 
with identical dipole moments /i we thus have (cf. Eq. 0): 

aim) , . 

Vc OC ^ . 21 



In Eq. ^ we have to insert the functional dependence /i(m) of the aggregate 
dipole moment on the aggregate mass. As shown in the previous section, 
this dependence has been found in numerical simulations ( [NuBOLD and 



Glassmeier, 2000| ; |iNUBOLD, to be fi{N) = ij,{m) oc m'^ with 7 = 0.63. 



In combination with Eq. ^we obtain oc mP'^^-m^^^"^ = m°'^^, i.e. u ~ 0.13. 

In contrast to diffusion-limited aggregation, the collisional cross section of 
magnetic grains is greatly enhanced over the geometric particle size. Again, 

^In the regime of free molecular flow, the mean free path of the molecules of the 
embedding medium (gas) is greater than the size of the dust particles under consideration. 
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for perfectly aligned dipoles we have: 

(Tc oc (m^)*^^ , 

= (m°-63)'/', (22) 

thence rj ^ 0.84. 

Since the foregoing derivations of u and t] imply maximum magnetic in- 
teraction, these values have to be considered as strict upper limits. To treat 
this problem in a more quantitative manner, we shall include the increas- 
ing dipole relaxation time T(m), i.e. the characteristic mass-dependent time 
scale for dipole alignment. 

The relaxation time r associated with a rotation of dipoles, e.g. through 
an angle = 7r/2, can be estimated from: 




where uj, /, and M denote the magnetic angular acceleration, the moment of 
inertia of the aggregate, and the magnetic inter-particle torque, respectively. 
Thus we obtain r oc V 7722-(i-0-63) _ jji^ -37 _ jf include this aspect in our 
calculation of the exponents u and t], we find the alignment-corrected values 
u* and f]* by replacing //^ in Eq. |^with ii^^^-^'^ = /i^'^^ to account for imperfect 
alignment. The results are: 

u"" ^ 0.01, (24) 
ri* ^ 0.68 . (25) 

Insertion of rj* and u* into Eq. ^ yields z ~ 3.2. This value for the 
dynamic exponent should characterize the kinetics of magnetic dust aggre- 
gation. It is the direct consequence of the combined effects of enhanced 
collisional cross sections and accretional remanence. 

4 Simulations using SAND: Two test cases 

A thorough analysis of the intricate microphysics of magnetic aggregation 
requires numerical tools which are able to treat both long-range magneto- 
static and short-range surface forces for large particle systems. To this end. 
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we devised a new coagulation code called SAND. Please refer to Sect, ^for 
a detailed description of the program which was also used for additional nu- 
merical simulation of magnetic coagulation. The results of this work will be 
the subject of a subsequent paper. In this section, we present two test cases 
which are meant to illustrate the abilities of SAND. 

4.1 Cross sections 

To test the validity of our theoretical treatment of magnetic aggregation, we 
performed detailed one-on-one scattering simulations, in which one dipole 
is magnetically deflected or attracted in the field of an identical, initially 
motionless, counterpart. 

For a detailed study of magnetic cross sections, we systematically varied 
the impact parameter around the value 6max as given by Eq. [TB|. We are 
thus able to determine the maximum value for collisions at different initial 
relative velocities v^o, which were taken from a Maxwellian distribution at 
100 K {0.3 -vth < <3-Vth). 

[Figure 2 about here.] 

[Figure 3 about here.] 

For all configurations [voo, b] we randomly chose the initial orientations 
{9i, (pi) of the two scattering dipoles, which were placed at a relative distance 
as given by Eq. |^. Figure ^ shows three example trajectories. Only the 
first scattering calculation with an impact parameter h = 0.95 ■ bmax leads 
to a "hard" collision of the grains. The second calculation with an impact 
parameter b = fomax does not lead to a "hard" collision. Instead we see that 
the two grains almost entirely exchange their energy and momentum. The 
impacting grain is moving only very slowly after the collision, while the grain 
initially at rest leaves the scene at a speed comparable to the impactor's initial 
velocity. In Fig. ^, different velocities are represented by different spacing 
of dots along the trajectory. For the third calculation with b = 1.05 ■ 6max, 
the grains do not hit each other and there is only partial transfer of energy. 
These three calculations therefore show good qualitative agreement between 
the predicted cross section (c.f. Eq. 0) and the numerical result. 

Figures @(a,b) analyze the results of a large number of simulations in a 
quantitative way, using color coding to account for orientation dependent 
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collision probabilities. While the detailed values are generally smaller than 
those obtained previously ( |Nubqld and Glassmeier, 1999| , pOOOj ), the 
qualitative agreement is good. At vth = ^/kBT/m we find: 

— {voo = Vth) ~ 130 . (26) 
a 

The discrepancy between the two results presented in Eq. |^ and Eq. ^ 
can be interpreted as a direct measure of the effect of the degree of magnetic 
alignment, although it should be noted that a small part of the discrepancy is 
due to the change in temperature (80 K vs. 100 K). According to Eq. |l^, am 
scales with T^^/^. The temperature-corrected value for am in Eq. |2^ should 
therefore be approximately 130 -O.S"^/^ ~ 151. The numerical values equally 
validate Eq. |T^. As had to be expected, the theoretically derived value for 
am - the dashed line in Fig. ^ - is too big, since only maximum interac- 
tion enters the central-force approximation. Orientation averaged collision 
probabilities reach 100% at smaller impact parameters that can be estimated 
from the outer radius of the white annuli, Rw This can be accounted for by 
introducing an empirical correction factor c = {Rw/RdashY in Eq. [1^: 

am = c-37i ( . 27 

By inspection of Fig. H we estimate c to lie between 0.8 and 0.9, i.e. the 



approximation error in Eq. 14 is of the order of 10-20%. We do not see a 



strong velocity dependence of this correction factor. 



[Figure 4 about here.] 



4.2 Restructuring and contact forces 

The long-range interactions between magnetized dust grains give rise to an 
increase in inter-particle collision velocities which have to be accounted for in 
realistic coagulation models. Figure ^ illustrates the importance of including 
configurational fiexibility and means for aggregate restructuring in coagula- 
tion models. Model calculations of dust coagulation usually dwell on very 
few assumptions, sometimes neglecting the aggregate structure altogether. 
The simplest approach treats grain growth in analogy to the coagulation of 
liquid drops, where the combination of two drops leads to the formation of 
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one - bigger - drop. Next in complexity are "hit- and- stick" models (e.g. 
Weidenschilling and Ruzmaikina, 1994| ; [Ossenkopf, 1993| ; [Meakin 



AND DONN, 1988D where grains are assumed to stick where they hit, whence 
the name, and do not influence the aggregate structure. While the particle 



trajectories may be either ballistic (e.g. [Ossenkopf, 199^ )) or stochasti- 
cal, i.e. pseudo-ballistic with random changes of direction on a diffusion- 
controlled time scale (|Kempf, 1997| ), there is no aggregate restructuring 



during collisions. Some models include sticking probabilities ([Ossenkopf 



1993| ) whereas others prescribe fixed values for the fractal dimension at dif- 



ferent stages of the simulation ( [Weidenschilling et al., 1988[ ). A realis- 
tic treatment of aggregates as solid bodies has been given by [Richardson 



The detailed analysis of inter-particle contact forces by POMINIK and[ 
I'lELENS (1997[ ) allows a more consistent approach. Clusters of particles 
are no longer treated as separate objects that must be parameterized with 
respect to their physical properties, rather they preserve their identity and 
are considered under the influence of all forces transmitted through each 
grain's contacts with adjacent particles. Especially when strong inter-particle 
forces are involved, grain-grain collisions can be violent and produce intricate 
results. This is shown in Fig. ^ where the linear aggregate is restructured to 
a T-shape after the fast impact of another grain. 



5 Conclusions 

The material presented in this paper lays the foundation for a detailed anal- 
ysis of the effects of magnetic aggregation for the accumulation of solid mat- 
ter in the early solar system. Generalizing previous work by [JNUBOLD AND 



Glassmeier (19"99[ , [2000| ), we have derived an analytical expression for the 



coUisional (capture) cross section of magnetic dust grains. As pairs of dipolar 
magnetic particles tend to line up in pole-to-pole configuration, aggregates 
of such particles "preserve" the magnetic moment of their constituents, thus 
displaying "accretional remanence" . The combination of enhanced cross sec- 
tion and accretional remanence is responsible for accelerated growth dynam- 
ics. This fact is quantitatively described by a dynamic exponent, z, such that 
the temporal dependence of the mean aggregate mass, m, takes on the form 
of a power law, i.e., m oct^. We have shown that the dynamic exponent takes 
on a value z ~ 3.2 for magnetic aggregation, whereas z = 2 for aggregation 
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processes with non-magnetic particles. 

Part of our theoretical approach has been validated numerically as a 
means for introducing a new coagulation code which we devised for the de- 
tailed analysis of magnetic coagulation. Its architecture allows the study of 
systems of particles interacting via long-range magnetostatic coupling as well 
as through short-range elastic forces. Aggregate restructuring is included. 
We provide evidence that such an approach - although computationally 
rather costly - is required for physically realistic modelling. A forthcoming 
paper will present extensive numerical studies of magnetic aggregation. 
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implement grain rotation, Josef Schiile for help with optimizing the code. 
CD acknowledges financial support from Pioneer grant 600-78-333. HN ac- 
knowledges financial support from the Studienstiftung des Deutschen Volkes 
and from DFG grant GL 142/13-1. 

A The SAND code 

SAND (Soft Aggregate Numerical Dynamics) is a code to compute the ag- 
gregation of dust particles in a detailed way, treating the microphysics of the 
contacts between dust particles, long range forces like magnetic dipole and 
electrostatic interactions, and external forces like external magnetic or grav- 
itational fields. This is a new combination which in this way has not been 
used before. The precursor of SAND was the 2D code used by POMINIK 
AND TiELENS (1997| ) to study the behavior of dust aggregates in collisions, 
but without long-range forces. 

SAND treats aggregates of dust grains as N-particle systems. Each in- 
dividual grain is free to move under the forces acting on it. Such forces 
can be due to external fields, to long-range grain-grain interactions, and to 
mechanical force transmission through physical contacts between grains. 

SAND currently treats only spherical particles, but an extension to sphe- 
roids or ellipsoids would be possible. 

SAND was original written to treat contact forces only, and therefore 
was best suited to study the collision of aggregates. Full aggregation cal- 
culations were not possible because the detailed numerical treatment of the 
contacts imposes very short time steps, of the order of 10~^^ seconds in or- 
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der to resolve the vibrational motion of grains in contact. Long-range forces 
considerably shorten the growth-time for an aggregating particle system. It 
therefore becomes feasible to conduct coagulation calculations with the code. 
On the other hand, implementing long-range forces also brings up the fol- 
lowing challenges. 

• For long-range forces, every particle in an A^-particle numerical simu- 
lation interacts with all other particles, thus leading to an N"^ scaling 
for the computation time. If only contact forces are involved, the code 
is linear in N. 

• Some kind of boundary condition has to be introduced in order to model 
the interaction of "edge particles" , i.e. particles close to the edge of the 
spatially limited simulation volume. For forces of virtually unlimited 
range, this can be a serious complication. 

• As interaction forces are strongly dependent on the inter-particle dis- 
tance, interaction between grains that are far apart will be treated in 
too much detail - numerically speaking - if only one global simulation 
time step is used. This problem is of particular importance in the case 
of SAND, since strong contact forces require very small time steps in 
the numerical treatment of the differential equations involved. 

SAND tackles these problems by means of a list-based interaction calcu- 
lation for long-range forces and a multiple time step algorithm, that separates 
short-range contact forces from long-range dipolar interaction. 

In the following we describe briefly the equations solved by SAND. Two 
test calculations are discussed in Sect. ^. More detailed simulations of mag- 
netic accretion will be the subject of a forthcoming paper. 

A.l Grain properties 

Each individual grain i is described by material properties (specific density 
Pi, surface energy 7^), radius Oj, mass rrii, moment of inertia Jj, magnetic 
moment vector /Xj. The state of each particle is given by a location r,, a 
velocity Vi, rotational state g^, and by a rotation velocity cjj. 
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A. 2 Linear Motion 



The differential equations describing the hnear motion of each grain are sim- 
ply 

r-i = Vi (28) 

'^j -^i, contact ~l~ -^i, magnetic • (2^) 

-Pi.contact is the sum of all linear elastic forces transmitted through contacts 
involving grain i. -Fj^magnetic is the sum of all magnetic interactions of grain i 
with other grains and an external magnetic field. In a similar way, gravitation 
and electrostatic forces could be added, but this is not relevant for our current 
study. 



A. 3 Rotation 



The rotational state is a four-vector which contains the rotation param- 
eters Qi^i, gi^2, Qi,3, and Qi^^ which can be used to describe the rotation of a 
solid body by differential equations, avoiding the singularities characteristic 
for Eulerian angles and similar 3-parameter representations. We follow the 
formalism by |WHlTTAKER (19041) . 

Grains initially enter the calculations with parameters Qi{to) = (0, 0, 0, 1). 
At a later state, g^it) = (f?i,i(t), 2(^)1 f?j,3(^), f?j,4(^))- The new orientation 
can be computed from the old orientation by a single rotation using the 
rotation matrix 



/ Qi,lQi,l — Qi,2Qi,2 — 



'^{Qi,iQi,2 + Qi,3QiA) 

'^{Qi,lQi,3 ~ Qi,2Qi,4:) 



'^{Qi,iQi,2 - Qi,3Qi,A) '^{Qi,iQi,3 + Qi,2Qi,A) 



\ 



— Qi,lQi,l + Qi,2Qi,2 — 
Qi,3Qi,3 + 

'^{.Qi,2Qi,3 + Qi,iQiA) 



'^{Qi,20i;3 ~ Qi,10i,'l) 



^Qi,lQi,l ~ Qi,2Qi,2 + 

Qi,3Qi,3 + Qi,iQi,A / 

(30) 



where we have omitted the time arguments for brevity. This matrix can be 
used to rotate any vector properties of the particle, such as the magnetic 
moment into the current state. 
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The differential equations for the rotation parameters are given by 

Qi,2 = {QiA^i,2 - ft,3^i,l + e'i,l^^j,3)/2 
Qi,3 = {Qi,2i^i,l - Qi,li^i,2 + Qi,4^i,3)/'2 

QiA = -{Qi,3^^i,3 + Qi.2^i.2 + Qi^\0Ji^\)l1 . (31) 

The rotation parameters also must fulfill a normalization condition, i.e. 

qIx + QI2 + QI3 + q\a = 1 • (32) 

Solving the differential Eq. |31| introduces numerical errors into the rotation 
parameters. After each successful timestep, we therefore re-normalize the 



rotation parameters in accordance with Eq. 32 . 

The differential equation for the change of the angular velocity u is simply 

— -^^i, contact ~l" -^^i, magnetic ; (^3) 

where -Mj^contact is the sum of all elastic torques transmitted through contacts 
involving grain i. Afj^magnetic is the torque due to all magnetic interactions 
of the grain i with other grains and with an external magnetic field. 

A. 4 Magnetic forces 

The basic equation for the magnetostatic interaction between two dipolar 
particles has been given in Sect. |3.1| (cf. Eq. |^). For notational simplicity, 
we let the constant factor {fiQ/in) = 10~^ be equal to 1 for the remainder of 
this section. The derivation of the force acting between two magnetic dipoles 
i,j is straightforward. In three-dimensional Cartesian coordinates, we have: 

F,^ = [fi^V) Bi, . (34) 

In Eq. |3^, the double subscript ij describes the influence of particle i on 
particle j, i.e. at the location of particle j. Individual particles are referred 
to by single subscripts, e.g. i. 

Carrying out the multiple differentiations in Eq. and introducing the 
relative distance vector, Vij, of length r we have: 



F = 1 



r 



2 



/x-, (35) 



with Ml = ■ Vij'^ + rij ■ + (/x^ ■ Vij) ■ I, (36) 
and M2 = ■ Tij'^ . (37) 
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In Eqs. ^ and the superscript T denotes a row vector. The matrix X 
in Eq. ^ is the identity matrix. All matrices contained in Eqs. p5|-|37| are 
symmetric, i.e., for actual computation only six out of nine elements have to 



be calculated. Trivially, F 



-F 



Apart from translational forces, dipoles will also exert torques Mi^t on 
their counterparts. Each contribution iVfj^ int, i.e. the torque acting on dipole 
j in the magnetic (dipole) field due to particle i at the location of j, is given 
by 



jj'jint 



fJ'i X Bi 



3^ 



Mi 



(38) 
(39) 



In addition to mutual interaction (cf. Eq. ^9]), magnetic grains can in- 
teract with an external magnetic field Bcxt- As a result of an undamped 
interaction the perpendicular part of the magnetic dipole moment vector of 
each grain, /z* = ^ij ^b-> ^ill oscillate around the direction of the external 
field. The magnetic torque Af j ^xt causing this oscillation is given by 



M 



j,ext 



\M 



j,ext I 



•"cxt 



(40) 
(41) 



In the scope of SAND, the external magnetic field will be constant in time 
and - without loss of generality - directed along the z-direction. 



A. 5 Contact forces 

When two grains collide, they can form a contact. Forces are transmitted 
through such contacts, and the magnitude of these forces largely determines 
the structure and mechanical properties of an aggregate. A contact between 
two grains has an equilibrium state in which the attractive forces are bal- 
anced by elastic forces due to compression of the material near the contact 
( IDOMINIK AND TiELENS, 1995| , |1996| , p,997| ). Dynamic forces on the grains 



lead to additional stresses in the contact area which are responsible for con- 
tact forces acting on the two grains. As long as the forces do not exceed 
the elastic regime, the contact forces can be calculated as a function of the 
displacement of the grains from the equilibrium contact position. Therefore, 
the problem of computing contact forces is reduced to computing these dis- 
placements. If the elastic limits of the contact are exceeded, the contact 
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breaks or moves, and energy is dissipated in this process. In the following 
we describe the procedure which is used in SAND to keep track of contact 
locations and displacements. 

When a new contact between grains i and j is formed, we define the 
following vectors describing the position of the contact on each grain. Let 
c denote the position where the grain surfaces touch. In the program, this 
defines a new contact structure which is kept in a list. We define a vector 
Pi = c — Vi which points from the center of grain i to the contact point 
and a vector P2 = c — Vj which points from the center of grain j to the 
contact point. We also define a unit vector Q which is perpendicular to the 
line connecting the two grain centers: |Q| = 1 and Q ■ {vi ~ Vj) = 0. 



At a later stage during the calculation, the two grains i and j have moved so 
that their positions and/or rotational states are different. 

If a contact were treated as being entirely rigid, the grains could only 
move such that the contact remains in its original position on the connecting 
line. This corresponds to a treatment of aggregates as rigid solid bodies 
( IRlCHARDSON, 1995|) . However, both grains are still treated as individual 
particles and they can move in a different way. If the contact is not lost, 
the material near the contact is strained in order to accommodate the small 
differences from rigid body motion. We can now use the vectors Pi, P2, 
and Q to compute how much the material has been strained. Figure |^ shows 
the definition of these vectors. The displacements can then be calculated in 
the following way. First we rotate the vectors Pi, P2, and Q into the new 
orientations of the grains, using the current rotation matrices 



Pi = 7^,Pl , P2 = T^jP2 , Qi = n^Q , 92 = n,Q . (42) 



The total linear displacement between the contact pointers is the vector 



We define a unit vector pointing from the center of grain j to the center of 
grain i 



[Figure 5 about here.] 



(43) 



(44) 



e± = 
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and compute the linear displacement in the direction of the inter-grain axis 
5± (pulling degree of freedom) and the displacement in the contact plane d\\ 
(sliding degree of freedom) with 

S± = {S- e^)e^ (45) 
Sii = S-S±. (46) 

Due to rolling of the surfaces over each other, the center of the contact area 
is displaced from the inter-grain axis by 

^ (Pi - P2) X jPi X P2) 

Finally we can compute the relative angular displacement around the inter- 
grain axis a\\ (twisting degree of freedom) and around an axis in the contact 
plane but perpendicular to the ^ vector: a± (rolling degree of freedom). 

a\\ = arcsin(qr;^ x qfg) ' (48) 
a± = aiccos{pi ■ (— P2)) • (49) 

From these displacements, the forces and torques acting on the grains 
can be directly computed using the formalism described by [DOMINIK and 
I'lELENS (1997| ), with the following modifications: 

• Heim et al. (1999) and [Blum and Wurm (2000| ) have shown that the 
critical rolling displacement ^crit is of the order of 10-20A instead of lA 
as used by pOMlNiK and Tielens (199^) . For the calculations in this 



paper, we have used ^crit = lOA. Rolling is the dominant restructuring 
process in the simulations in this paper, and a proper treatment is 
important here. The measurements by Blum and Wurm have been 
restricted to Quartz particles, so the application to metallic or ceramic 
particles introduces some uncertainty. 

I — I 

• PoPPE et al. (2000) have shown experimentally that the critical ve- 
locity for sticking can be considerably higher than the value predicted 
by [Chokshi et al. (1993|) . This has to be interpreted as an additional 
channel of energy dissipation which is not yet well understood. Since 
there is no good theoretical description available, we have still used 
the lower energies and velocities for sticking. In the current study this 
effect is probably not of critical importance since the attraction be- 
tween grains is dominated by magnetic forces and generally leads to 
immediate sticking. 



25 



A. 6 Restructuring: making, breaking, moving contacts 



The algorithm described in the previous section is the one used while contacts 
are in the elastic regime. In this regime, the aggregate keeps its basic shape 
except for minor bending and streching. When the dynamic forces acting on 
grains in contact become too big, the contact will start to move over the two 
grains by rolling or sliding, or it will break. Within each time step, SAND 
assumes all contacts to be stable and elastic, in order to make it possible 
for the integrator to work on a set of continuous functions. This assumption 
limits the time step possible for SAND to the period of the fastest contact 
vibration. After each timestep, SAND checks if any of the contacts are 
strained too much. If this is the case and the grains are too far apart to 
stay in contact, the contact is broken and removed from the internal list. 
If the grains are still close enough to touch, the contact is moved to a new 
equilibrium location, releasing the stresses that prompted the motion. Then 
new values for the contact vectors Pi, P2, and Q are computed and the time 
integration continues. Also after each time step, the integration volume is 
checked for grains which have approached each other so that a new contact 
has to be establised. 

A. 7 List algorithm 

An A^-particle code with long-range forces can only be numerically efficient, 
if we can cut down the O (A^^)-dependence to a more acceptable 0{N)- 
dependence. To within a numerical factor, this obviously represents the 
lower bound for the amount of work required to process all particles. 

In Eq. 1^, we defined a thermally induced interaction limit for dipolar 
magnetic coupling. We can thus effectively reduce the long-range interaction 
of a given dipole j to that with a restricted set of neighbor particles. The 
computational problem is now reduced to implementing an efficient algo- 
rithm for list construction and maintenance. We have used a cell subdivision 
( |Rapaport, 1995]) adapted for the use with magnetically interacting grains. 

We use Eq. |^ as a definition for the cut-off distance, r^, for our long-range 
forces. The simulation space is subdivided into cells whose (identical) edges 
all exceed Vc in length. Interaction is only possible if two grains are either in 
the same cell or in immediately adjacent (nearest neighbors) cells. In three 
dimensions, each cell has 26 neighboring cells. Due to reasons of symmetry 
only half of these need be considered. Note that this method only works 
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properly if the size simulation volume exceeds Ar^. Since we do not know in 
advance how many particles will occupy a given cell at any instant, linked 
lists ( [Rapapqrt, 19951) are used to associate particles with the cell they 
reside in. 

The cell method can be considerably refined by a neighbor-list method. 
For cells with an edge length , on the average only a small fraction of 
the examined particles lie within interaction range. In order by make use 
of this reduced neighborhood size, we have to construct a list of interacting 
pairs from the set of particles found by the cell method. To this end, we 
replace the cut-off distance Vc by Vc + Ar with Ar > in order to allow 
the neighbor-list to be useful during several numerical time steps, Atj. The 
concise value of Ar is related to the rate at which the list has to be updated 
in an inverse fashion. The list is updated whenever the cumulative maximum 
velocity i^max exceeds a certain limit: 

v^, , Ar , , 

2^ l^'^.ma.l > ■ (50) 

This is a rather conservative criterion since it sums contributions from dif- 
ferent particles. In practice, Ar-values between CSr^ and 0.4rc have proved 
useful. For the evaluation according to Eq. ^ only unbound grains are con- 
sidered. Bound particles are subject to strong elastic forces ( [DOMINIK and 
I'lELENS, 1997] ) that cause rapid oscillation. 

In order to include long-range forces on particles close to the edge of 
the simulation volume, "wraparound" boundary conditions have been imple- 
mented into SAND. 

For an ensemble of 100 magnetic particles, the combined cell-based neigh- 
bor-list method leads to a decrease in computation time by a factor of 30 to 
40. Since our approach is validated not only on a numerical basis but also 
with respect to physical considerations we have implemented an efficient way 
of dealing with magnetic forces in an A^-particle simulation. 

A. 8 Multiple time scales 

Based on the considerations in the last section, dipolar magnetic interaction 
of thermalized grains, i.e. grains that are subject to random thermal motion 
and rotation, can be viewed as moderately long-ranged with cut-off distance 
Tc. On these grounds the complete set of A^(A^ — l)/2 pair interactions can 
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be divided into a so-called cardinal set which contains all pairs for 
which dij < Tc, and the so-called inessential set for which interactions are 



neglected ( |Streett et ai, 1977| ). The population of both the cardinal and 



the inessential sets change with time. 

The multiple time step (MTS) approach focuses on the point that the 
cardinal set can be subdivided into classes of primary and secondary neigh- 
bors based on a distance criterion for potentials with strong r-dependence. 
It is the rapidly changing force due to the primary neighbors that sets an 
upper limit for the time-step At in the simulation. In the MTS method, 
only the contributions of these primary neighbors are updated at every step. 
Longer time-steps are used to calculate the evolution of the secondary force, 
i.e. the contribution due to the secondary neighbor particles. To achieve 
this, the secondary force at time t = to + kAt is calculated from a Taylor 
series expansion based on information at a given time to: 

F^ito + kAt) = FM + Vf^H^o) ^ (51) 

ml 

m 

where the superscript (m) denotes the differentiation with respect to time. 
A similar equation could be written down for the torque, M. Equation |51| is 
referred to as the mth-order Taylor series. The maximum value fcmax for the 
parameter k has to be chosen in such a way that the gain in computational 
efficiency balances the loss in numerical accuracy. During the simulation, k 
runs from 1 to /cmax- At that point a new Taylor series is calculated according 
to Eq. |5l| and that step becomes the new to- 

While using additional memory for the storage of Taylor expansions, this 
method can achieve an additional increase in computational speed for all 
but the smallest systems up to a factor of five ( [Streett et al., 1977| ). The 
authors also report on the numerical error as a function of the expansion 
order m and the parameter fcmax- 

In our coagulation model, an MTS approach is already inherent in the 
numerical setup. In the scope of SAND, we identify the elasto-mechanical 
contact forces with the primary force and the magnetic forces with the sec- 
ondary force. Since we use an adaptive step size algorithm to integrate the 
equations of motion, the step length At is automatically adapted to the 
"physical" situation based on a prescribed absolute error criterion. As soon 
as the first inter-particle contact is formed, the average step size drops by a 
factor of 10^ due to the strong elastic forces acting on only a few particles. 



28 



At this point we call upon the MTS method, which has been implemented 
as an optional feature into SAND. 

Although straightforward mathematically, the calculation of the Taylor 
expansion for magnetic dipole-dipole interaction is a rather cumbersome task 
for Fij = Fij{iJ,^, fj,j, Vij). As before, the vector Vij denotes the relative posi- 
tion of the two dipoles involved. In what follows we shall drop the superscript 
{ij) for simplicity's sake. Unless otherwise noted all vector quantities are as- 
sumed to originate from particle i and to be directed toward or acting upon 
particle j. The basic formula involved is: 

^ " dr dt + dfi, dt ' ^^^^ 



k=i,j 



which is the total time derivative of the dipolar force. Since SAND solves 
the full set of equations for particle motion including rotation, the vectors (Xf. 
are simply given by ■ cj^ if ujk denotes the angular velocity of the grain. 
V is the translational velocity of the particle. The matrices Jr and Jk,^ are 
the Jacobian matrices. 

The Taylor expansion for the magnetic torque has to incorporate both 
internal and external contributions. As before, we drop the subscript ij. 
The expression for the first time derivative of the torque Adj due to internal 
interaction then reads: 

M,,.t = ^ X S + X ^ . (53) 

The expression for the dipolar magnetic field B has been stated before (cf. 
Eq. In analogy with Eq. ^ the time variation of the magnetic field, 
dB / dt, can be determined from 

B = Jrr + Ji^^ /X, . (54) 

The external contribution for each magnetic grain is simply given by: 

Mj^ext = IMjI • X ^ext ■ (55) 

We therefore use the following expression for the total time derivative of the 
magnetic torque on a given particle j which can be derived by combining 
Eqs. 

Mj = fXjX B + X {J, r + Ji^^ fii) + l/Xjl ■ c^j x Bcxt • (56) 
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In our numerical work, Taylor series expansions have been constrained to first 
order due to the rather complicated structure of the derivatives. These are 



recalculated at to according to Eqs. B2 and 56. During the following fcmax time 



steps, which are governed by the high numerical demands of strong contact 
forces, all magnetic forces are calculated using Eq. ^ and its counterpart for 
the torques. In practice, the choice fcmax = 70 has yielded satisfactory results 
both in terms of computing time and physical accuracy. 
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FIG. 1. Accretional remanence: The increase of cluster magnetic mo- 
ment with the number of individual magnetic grains in numerical simulations 
(INUBOLD AND Glassmeier, 1999| , |2000|) . The curve illustrates the rela- 
tionship fi oc N^-^^^ which shows the strong tendency of magnetized grains to 
aquire a state of minimum magnetostatic energy (cf. Eq. i.e. a pole-to- 
pole configuration for each pair of grains, and the misalignment effect caused 
by thermal agitation as well as the microphysics of the coagulation process. 
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FIG. 2. Three example trajectories of numerical collision experiments 
with SAND. The particles are 20 nm iron grains with a magnetic moment of 
200Am^kg~^ and thermal rotations derived from a Maxwellian distribution 
at 100 K. Initial velocity is Uth/16. The impact parameters are 0.95, 1.0, 
and 1.025 times the analytically derived maximum impact parameter femax, 
respectively. Note that the size of the dots does not indicate the size of the 
grains which would be hardly visible on this scale since fomax ~ 250i?. 
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(a) Low velocity part. 



(b) High velocity part. 



FIG. 3. Numerical validation of Eq. |T4[ Collision probabilities for dif- 
ferent relative velocities and impact parameters are orientation averaged and 
color-coded. Black means that for no relative orientations of dipoles a colli- 
sion could be detected. White means that all randomly chosen orientations 
resulted in a collision. For one specific velocity, all the area outside the black 
annular region should be black as well, the same holds for the area inside 
the white annular region, which should be all white. Please note that the 
color bar has been compressed for the high velocity regime, with black coding 
meaning p < 60%. The dashed lines denote the theoretically derived cross 
section according to Eq. 



38 



i 

I 


1 

1 


! 


• - 

I 


L 1 . 


s 

8 


■t 



FIG. 4. Example of a collision between a chain and a single particle 
which illustrates that the detailed treatment of contact forces is important to 
determine the correct structure of an aggregate. The impact of the incoming 
grain moves the grain at the end of the chain. The resulting T-like structure 
would not have resulted from a simple hit-and-stick approximation. 
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FIG. 5. Geometrical treatment of contact displacements. Left panel: 
Contact between two grains in equilibrium (top) and in stressed configuration 
(bottom). Right panel: Blow-up of the contact region with definitions of the 
displacement vectors. See text. 
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